pacman::p_load(dplyr, purrr, ggplot2, GGally, tidyr, sf,
factoextra, ggdendro, cluster,plotly, paletteer,
patchwork, umap, MASS)
options(readr.show_col_types = FALSE)
set.seed(42)
theme_set(theme_bw())
create_graphs = TRUETarea 1: Cluster Espaciales
Métodos de aprendizaje de máquinas en Data Science
Enunciado General
El territorio de Chile esta subdividido en unidades administrativas, como regiones, provincias y comunas. La definición de estas unidades actuales obedece a diversos motivos históricos, políticos y circunstanciales, por lo que no existe homogeneidad entre ellas y existen comunas de diversos tamaños, poblaciones y condiciones de vida. La Subsecretaría de Desarrollo Regional y Administrativo de Chile ha decidido que para optimizar el diseño y ejecución de políticas públicas necesitan crear unidades administrativas que tengan un mayor grado de homogeneidad interna, también conocido como cohesión, y que estas regiones estén cohesionadas también en el espacio.
Para estos efectos, la SUBDERE lo ha contratado a usted como consultor para que diseñe un algoritmo de clusterización que permita crear regiones en el territorio que cumplan con estas características. Para cumplir esto, se le otorgará una base de datos de todas las manzanas de la comuna de Las Condes (~1.600), donde cada manzana contiene 39 atributos que reflejan diferentes dimensiones demográficas y de condiciones de vida, según lo reportado en el censo 2017 y lo calculado por la plataforma de Bienestar Humano Territorial.
El objetivo principal de este encargo es crear un programa computacional que permita crear clústers de manzanas que sean similares entre sí para reemplazar a los 13 cuadrantes actuales.
Análisis Descriptivo de Datos
A continuación se procede a realizar un proceso de análisis desciptivo general de los datos a modo de entender su distribución y relaciones entre ellas.
Lectura de Datos
df <- read.csv("data/datos_t1_centroides.csv") %>%
dplyr::select(-1)%>%
janitor::clean_names()
# glimpse(df)La base de datos llamada df tiene 1661 registros correspondientes a centroides de manzanas censales de la comuan de Condes, cada uno de estas manzanas manatiene 40 atributos (columnas) que corresponden variables de identificación, valores de indicadores de Bienestar Territorial (IBT), características sociodemográficas y coordenadas espaciales (UTM 19S), como a continuación se detalla:
- Variables de Identificación: cod_com, comuna, id_manz, manz_en, zona
- Valores de Indicadores de Bienestar Territorial: ibt, dim_acc, iav, icul, idep, isal, iser, ise, dim_amb, iata, icv, dim_soc, ivi, isv, iej, irh, iem, ipj, dim_seg, igpe, igpr, ilpe, ilpr
- Variables de sociodemográficas: area, total_v, hog_n, personas, e0a5, e6a14, e4a18, e0a18, e15a24, e65ymas
- Coordenadas : x, y
Para visualizar espacialmente la distribución estas manzanas se procede a generar un mapa, cuyo color de cada centroide (punto) corresponde al valor del indicador de Bienestar Territorial (ibt).
# trasformación de dataframe a Simple Feature (objeto espacial en R)
sf_ibt <- df %>%
st_as_sf(coords = c("x", "y"), crs = 32719) %>%
mutate(zona_m = substr(zona, 9, nchar(zona)))
mapview(sf_ibt, zcol="ibt", cex = 3)Distribución Variables de Identificación
Como se mencionó las variables de identificación (cod_com, comuna, id_manz, manz_en, zona) tiene como fin identificar cada una de las manzanas en un contexto de otra uninidad adminustrativa mayor, como lo es Zona Censal (zona) y comuna a través de su código (cod_com). De estas las únicas presentar variabiliad y aportarn información es la Zona Censal y la identificación de manzana (id_manz). Se procede a visualizar un gráfico que representa la cantidad de manzanas por zona censal.
Indicador de Bienestar Territorial
Distribución Variables Sociodemográficas
Personas por Manzana
Viviendas por Manzana
Hogares por Manzana
Población de 0 a 5 años
Población de 6 a 14 años
Población de 4 a 18 años
Población de 0 a 18 años por zona
Población de 15 a 24 años por zona
Población de 65 y más años por zona
Distribución indicadores
Indicador de Bienestar Territorial
Dimensiones
Indicadores de Accesibilidad
Indicador de Áreas Verdes (Aparte la escala difiere enormente a los demás dimensiones)
Indicadores de Accesibilidad
Indicadores Ambientales
Indicadores Socioeconómicos
Escolaridad Jefe de Hogar (se hace aparte, mantiene otra escala)
Limpieza de Datos y Selección
En esta sección se trabaja sobre el espacio de features de los datos, revisando sus distribuciones y a partir de esto tomando decisiones metodológicas respecto del tratamiento de datos perdidos, posibles transformaciones de variables y reducción de dimensionalidad.
En esta sección revisamos la distribución univariada de cada variable de indicadores y algunas relaciones entre estas, se llevan a cabo transformaciones de variables.
# tipo de variables (identificador, valores y caracteristicas)
cols <- colnames(df)
id_cols <- cols[1:5]
val_cols <- cols[6:28]
car_cols <- cols[29:38]Manejo de Datos Perdidos
En la base podemos encontrar información faltante en ciertas variables, en particular:
map_df(df, ~ sum(is.na(.x))) %>% keep(~.x>0)# A tibble: 1 × 11
iav icul idep isal iser ise ivi isv iej irh iem
<int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
1 180 180 180 180 180 208 166 166 180 180 180Las variables que presentan información perdida se relacionan a indicadores que calculan accesibilidad a servicios para la población o información relacionada a viviendas y personas. Por lo que al revisar cómo se comportan estas observaciones tenemos bajos niveles de población, muy cercano a cero:
NA_rows <- rowSums(is.na(df))>0
my_summarise(filter(df, NA_rows), car_cols)# A tibble: 10 × 9
var min q25 mean median q75 max NAs N
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 area 249. 2478. 20220. 4673. 10495. 411415. 0 208
2 e0a18 0 0 0.154 0 0 5 0 208
3 e0a5 0 0 0.125 0 0 4 0 208
4 e15a24 0 0 0.490 0 0 14 0 208
5 e4a18 0 0 0 0 0 0 0 208
6 e65ymas 0 0 0.837 0 0 12 0 208
7 e6a14 0 0 0 0 0 0 0 208
8 hog_n 0 0 1.33 0 0 43 0 208
9 personas 0 0 3.54 0 0 77 0 208
10 total_v 0 0 1.38 0 0 34 0 208Es por esto que podemos descartar estas observaciones del análisis. La razón de esta decisión viene del objetivo de análisis, de encontrar unidades homogéneas en el espacio de variables y cohesionadas en el espacio geográfico. Por ahora podemos omitirlas del análisis, sin ensuciar el espacio de variables, y luego de detectar los clusters homogéneos pueden asignarse mediante cercanía a estos grupos.
# omitir NA del analisis
df_na <- na.omit(df)Transformación de Variables
Al revisar la distribución de los indicadores (Section 2.4) encontramos que algunas variables tienen patrones evidentemente sesgados a valores altos o bajos (isv, o ivi vs icv) además de que existen diferencias en las unidades de medida, las cuales son todas positivas. A partir de esto se aplicarán transforamciones logarítmicas buscando suavizar las distribuciones y asemejarlas más a una campana normal.
Así mismo las variables de caracterización muestran un sesgo a valores bajos, por lo que también se le aplicarán transforamciones logarítmicas.
Aplicamos una transformación logarítmica a las variables más sesgadas:
log_trans <- function(x, ep = 1) log(x+ep)
log_vars <- c("icv", "idep", "iej", "iem", "igpe", "ilpr", "ipj",
"isal", "ise", "isv", "ivi", "iser")
log_vars <- c(log_vars, car_cols)
df_trans <- df_na
df_trans <- df_trans %>%
mutate(across(all_of(log_vars), log_trans))Distribución indicadores Transformadas
Indicador de Bienestar Territorial
Dimensiones
Indicadores de Accesibilidad
Indicadore de Áreas Verdes (Aparte ya que esta en otra escala)
Indicadores de Accesibilidad
Indicadores Ambientales
Indicadores Socioeconómicos
Escolaridad Jefe de Hogar (se hace aparte, mantiene otra escala)
Estandarización y Outliers
Luego, aplicaremos una estandarización para que los valores estén centrados en el promedio y que cada unidad represente una desviación estándar:
# seleccion de variables
df_coords <- df_trans %>% dplyr::select(x, y)
df_vars <- df_trans %>% dplyr::select(!any_of(c(id_cols, "x", "y")))
# std
preproc <- caret::preProcess(df_vars)
X_std <- predict(preproc, df_vars)Al revisar la distribución expresada en términos de desviaciones estándar podemos apreciar claramente valores más de 4 desviaciones estandar del promedio:
my_summarise(X_std, colnames(X_std)) %>% print(n=nrow(.))# A tibble: 33 × 9
var min q25 mean median q75 max NAs N
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 area -3.31 -0.630 3.37e-16 -0.152 0.474 4.25 0 1453
2 dim_acc -4.65 -0.484 -5.83e-16 0.0777 0.623 2.41 0 1453
3 dim_amb -1.66 -0.796 9.07e-17 -0.104 0.808 2.07 0 1453
4 dim_seg -5.15 -0.636 2.22e-15 0.0689 0.490 1.81 0 1453
5 dim_soc -5.75 -0.290 2.66e-16 0.227 0.592 2.04 0 1453
6 e0a18 -3.16 -0.679 4.05e-16 -0.120 0.619 3.15 0 1453
7 e0a5 -1.94 -0.614 -3.48e-16 -0.0803 0.581 3.29 0 1453
8 e15a24 -2.94 -0.689 5.78e-17 -0.101 0.608 3.98 0 1453
9 e4a18 -2.33 -0.669 -1.49e-16 -0.120 0.628 3.13 0 1453
10 e65ymas -2.74 -0.585 -1.50e-16 -0.0880 0.561 3.43 0 1453
11 e6a14 -2.62 -0.694 5.22e-16 -0.144 0.601 3.18 0 1453
12 hog_n -2.27 -0.696 -1.05e-16 -0.274 0.572 3.43 0 1453
13 iata -4.21 -0.494 -3.18e-16 0.0625 0.619 2.67 0 1453
14 iav -3.56 -0.457 2.11e-16 0.0557 0.561 4.13 0 1453
15 ibt -4.98 -0.328 5.51e-16 0.148 0.642 1.94 0 1453
16 icul -2.26 -0.812 -1.17e-16 -0.132 0.739 2.85 0 1453
17 icv -1.45 -1.16 8.65e-16 0.236 0.841 1.57 0 1453
18 idep -2.08 -0.756 -7.36e-17 -0.0233 0.436 3.97 0 1453
19 iej -5.00 -0.00247 1.45e-15 0.288 0.487 1.41 0 1453
20 iem -7.95 -0.437 -1.10e-15 0.0581 0.535 1.33 0 1453
21 igpe -4.87 -0.0447 -1.53e-15 0.303 0.477 0.777 0 1453
22 igpr -4.47 -0.578 -1.13e-15 0.0548 0.646 1.73 0 1453
23 ilpe -4.01 -0.386 2.13e-16 0.158 0.504 1.52 0 1453
24 ilpr -12.3 -0.0474 3.29e-15 0.213 0.488 0.818 0 1453
25 ipj -15.0 -0.341 -1.15e-15 0.201 0.616 0.616 0 1453
26 irh -6.07 -0.454 -4.81e-16 0.112 0.555 1.76 0 1453
27 isal -4.43 -0.708 -2.87e-17 0.0489 0.842 1.89 0 1453
28 ise -3.51 -0.674 -4.13e-17 0.0654 0.621 3.57 0 1453
29 iser -1.70 -0.693 1.14e-16 -0.178 0.563 2.74 0 1453
30 isv -8.70 0.0688 8.41e-16 0.381 0.381 0.381 0 1453
31 ivi -12.7 0.267 2.34e-15 0.267 0.267 0.267 0 1453
32 personas -2.63 -0.662 -5.29e-17 -0.239 0.591 3.45 0 1453
33 total_v -2.75 -0.676 4.74e-18 -0.285 0.591 3.38 0 1453Estos pueden ser descartados para evitar afectar la estructura de relación entre las variables:
# descartar outliers
out_rows <- (apply(X_std, 1, function(x) any(abs(x)>4)))
X <- X_std %>% filter(!out_rows)
coords <- df_coords %>% filter(!out_rows)Análisis de Correlación
A partir de esto podemos revisar las asociaciones entre los indicadores trabajados:
corr <- cor(X)# improved correlation matrix
library(corrplot)
corrplot(corr,
method = "square",
# type = "upper", # show only upper side
tl.cex =0.9,
number.cex =0.4,
tl.col = "gray50"
)
Se aprecian grupos de indicadores positivamente relacionados como las caracterizaciones de personas y hogares, asi tambien indicadores como dim_acc iser, iav, iata e isal. Las asociaones dnegativas ocurren principalmente entre indicadores de seguridad y los indicadores de accesibilidad.
Feature Extraction
En esta sección reducimos la dimensionalidad del espacio de features para facilitar la detección de clusters, con el objetivo general de reducir tiempo procesamiento y memoria;facilitar la visualización y eliminar atributos irrelevantes eliminando ruido.
Usaremos distintas técnicas las cuales serán comparadas en su capacidad de generar una proyección de los datos donde se detecten más fácilmente clusters.
PCA
Comenzamos por análisis de componentes principales:
# pca
PCA <- princomp(X, cor=TRUE)
features_pca <-
predict(PCA) %>%
as_tibble() %>%
janitor::clean_names() %>%
dplyr::select(1:2)Al revisar la varianza explicada por cada combinación lineal de variables tenemos que con los tres primeros componentes no alcanzamos a expresar el 40% de la varianza.
barplot(cumsum(PCA$sdev/sum(PCA$sdev)))
UMAP
Luego aplicamos UMAP:
data.umap <- umap(X)
features_umap <- data.umap$layout %>% as.data.frame()Escalamiento muldimensional
# Multidimensional scaling
d <- dist(X) # distancias euclidianas entre entidades
data.MDS <- cmdscale(d, eig=TRUE, k = 2) # k es el numero de dimensiones de salida
features_MDS <-
data.MDS$points %>%
as.data.frame() %>%
as_tibble()Escalamiento multidimensional no paramétrico:
# nonparametric MDS
data.nMDS <- isoMDS(d, k=2) initial value 25.560808
iter 5 value 20.530199
final value 18.912842
convergedfeatures_nMDS <-
data.nMDS$points %>%
as.data.frame() %>%
as_tibble()Comparación
Comparamos sus resultados en su capacidad de generar clusters mediante el estadístico de Hopkins, donde se necesita probar la hipótesis de la existencia de patrones en los datos contra un conjunto de datos distribuidos uniformemente (distribución homogénea). A continuación se visualiza su evaluación:
h_umap <- get_clust_tendency(features_umap, n = 30, graph = FALSE)$hopkins_stat
h_mds <- get_clust_tendency(features_MDS, n = 30, graph = FALSE)$hopkins_stat
h_nmds <- get_clust_tendency(features_nMDS, n = 30, graph = FALSE)$hopkins_stat
h_pca <-get_clust_tendency(features_pca, n = 30, graph = FALSE)$hopkins_stat
list(h_umap = h_umap, h_mds = h_mds, h_nmds = h_nmds, h_pca = h_pca)$h_umap
[1] 0.9701799
$h_mds
[1] 0.72542
$h_nmds
[1] 0.7969982
$h_pca
[1] 0.7566609Se precia como UMAP genera la mejor extracción de features para identificar clusters en el espacio de variables. Se puede apreciar gráficamente:
ggplot(features_umap, aes(V1, V2)) +
geom_point()
Problema 1: Cluster en espacio de variables
- Enunciado:
- Aplique un algoritmo de clustering y seleccione 13 clusters sin considerar las variables de posición de las manzanas. Muestre el resultado del proceso de clustering (grafique los puntos del espacio con sus respectivas etiquetas de clusters), y vea cuales son sus ventajas y desventajas para la creación de las manzanas.
Usando la proyección de datos generada por UMAP identificaremos clusters en el espacio de variables y revisamos su expresión espacial.
Comenzamos por Cluster Jerárquico y K-medias, para ambos calculamos primeramente la distancia eucledian
# calculamos la distancia euclideana
d0 <- dist(features_umap)Cluster Jerárquico:
# hacemos un modelo jerarquico con distancia completa
model_complete <- hclust(d0)
# obtenemos una sintesis del modelo
hgroups <- cutree(model_complete, k = 13)kmeans
# kmeans
kmod <- kmeans(features_umap, 13)
kgroups <- kmod$cluster
Discusión y Evaluación
Se aprecia como las formas se sobrelapan en ambos métodos en la mayoria de asignación de, siendo sin embargo el cluster jerárquico más limpio en los grupos identificados de acuerdo a una inspección visual.
Problema 2: Cluster en espacio geográfico
- Enunciado:
- Aplique un algoritmo de clustering y seleccione 13 clusters considerando solo las variables de posición de las manzanas. Muestre el resultado del proceso de clustering (grafique los puntos del espacio con sus respectivas etiquetas de clusters), y vea cuales son sus ventajas y desventajas para la creación de las manzanas.
Ahora repetimos el ejercicio solo usando el espacio geográfico de las manzanas, usaremos tanto Cluster Jerárquico como K-medias y los comparemos.
Primeramente se calcula la distancias eclediana por corrdenadas, esto nos da una relación de proximidad real
d1 <- dist(coords)Cluster Jerárquico:
hgeo <- hclust(d1)
hgeo_groups <- cutree(hgeo, k = 13)kmeans
kgeo <- kmeans(coords, 13)
kgeo_groups <- kgeo$clusterGráficamente tenemos:
Discusión y Evaluación
Se aprecian clusters más claros espacialmente, además de la diferencia entre grupos encontrados en ambos métodos. Nuevamente a través del proceso de inspección visual encontramos que existen patrones concentración espacial con mayor definicón en cluster jerárquico denominado hgeo_groups.
# trasformación de dataframe a Simple Feature (objeto espacial en R)
sf_clusters <- df_coords %>% filter(!out_rows) %>%
st_as_sf(coords = c("x", "y"), crs = 32719) %>%
mutate(kgeo_groups = factor(kgeo_groups),
hgeo_groups = factor(hgeo_groups))
mapview(sf_clusters, zcol="hgeo_groups", cex = 3,
col.regions=paletteer_c("grDevices::rainbow", 13),
layer.name = "hgeo_groups"
)Problema 3: Cluster en espacio de variables y geográfico
- Enunciado:
- Aplique un algoritmo de clustering y seleccione 13 clusters considerando las variables descriptivas y las variables de posición. Sin embargo, para este problema, se le pide que pueda integrar la distancia entre las manzanas como algo restrictivo. Específicamente, que la distancia entre dos puntos que pertenecen a un mismo cluster no sea mayor a 2,000 metros (en la medida de lo posible). Muestre el resultado del proceso de clustering (grafique los puntos del espacio con sus respectivas etiquetas de clusters), y vea cuales son sus ventajas y desventajas para la creación de las manzanas.
Finalmente usaremos un método de cluster jerárquico geográfico implementado por la librería ClustGeo el cual usa dos matrices de distancias, una en el espacio de variables y otra en el espacio geográfico. Estas son combinadas de forma convexa y usando el parámetro alpha que le da la importancia relativa a cada una de las distancias. Cuando este es 0 la distancia geográfica no tiene imporancia y cuando es 1 tiene toda la importancia.
Para elegir el valor de este parámetro alpha se calculan las homogeneidades relativas en cada matriz de distancia en los grupos encontrados, la idea es elegir un valor que no implique un sacrificio importante en ambos frentes.
Comenzamos por implementar el método usando las distancias geográficas crudas y luego lo implementamos usando una restricción de vecindad de menos de 2000 metros:
library(ClustGeo)
range.alpha <- seq(0,1,0.1)
K <- 13
cr <- choicealpha(d0, d1, range.alpha,
K, graph = FALSE)
cr$Q # proportion of explained inertia Q0 Q1
alpha=0 0.9737437 0.8094104
alpha=0.1 0.9719445 0.8207395
alpha=0.2 0.9708055 0.8410808
alpha=0.3 0.9699829 0.8507461
alpha=0.4 0.9646048 0.8743724
alpha=0.5 0.9626346 0.8829204
alpha=0.6 0.9637729 0.8816922
alpha=0.7 0.9441026 0.8986122
alpha=0.8 0.9295569 0.9070623
alpha=0.9 0.8963452 0.9197628
alpha=1 0.7512889 0.9300041plot(cr)
A partir de esto se define un alpha de 0.9 que le da mayor relevancia al espacio geográfico y no implica una pérdida importante en el espacio de variables:
# alpha que compensa perdida en homogeneidad en features (D0)
# con ganancia en homogeneidad espacial (D1)
alpha_geo = 0.9
tree <- hclustgeo(d0,d1,alpha = alpha_geo)
geoclust <- cutree(tree,13)Para el caso de restricción de vecindad repetimos el proceso:
library(sf)
library(spdep)
pts <- st_as_sf(coords, coords = c("x","y"), crs = 3857)
nb <- dnearneigh(pts, d1 = 0, d2 = 2000)
A <- nb2mat(nb,style="B")
diag(A) <- 1
D1 <- as.dist(1-A)
range.alpha <- seq(0,1,0.1)
K <- 13
cr <- choicealpha(d0, D1, range.alpha,
K, graph=FALSE)
plot(cr) # alpha buen compromiso0.5
Se escoge un valor de 0.5 que no genera mayor pérdida en ambos espacios.
alpha_nb = 0.5
tree_nb <- hclustgeo(d0,D1,alpha = alpha_nb)
geoclust_nb <- cutree(tree_nb,13)Los resultados se presentan a continuación:
Discusión y Evaluación
En este caso a través de un proceso de inspección visual de resultados espacializados, encontramos que cluster denominado geoclust que utiliza las distancias crudas, con un factor de relevancia de ellas de 0.9 tiene resultados más consistentes o definición de los patrones de concentración, que el método que utiliza una restricción de distancia (2000 metros). A continuación se visualiza el cluster elegido.